#include<iomanip>
#include"PMesh.h"
#include <iostream>
#include <fstream>

using std::ifstream;
using std::ios;



PMesh::PMesh(){
	elements=NULL;
	nodes=NULL;
	numOfE=numOfN=hsize=vsize=0;
	beta=dt=0;
}

PMesh::~PMesh(){}

PMesh::PMesh(const int &hsize,const int &vsize){
   this->hsize=hsize;
   this->vsize=vsize;
   this->numOfN = (hsize * vsize);
   this->nodes = new Node[numOfN + 1];
   this->numOfE = 2 * (hsize - 1) * (vsize - 1);
   this->elements = new PElement[numOfE + 1];
   beta=dt=0;
}

Matrix PMesh::assembleGlobalCond(const double &rho, const double &cp)
{
	Matrix g=Matrix::zeros(numOfN,numOfN);
	for(int i=1;i<=numOfE;i++){
		elements[i].constructCond(rho,cp);
	}
	
	int p=0;
	int r=0;
	
	for(int e=1;e<=numOfE;e++){
		for(int i=1;i<=numOfN;i++)
		{
			p=elements[e].indexOf(nodes[i]);
			for(int j=1;j<=numOfN;j++)
			{
				r=elements[e].indexOf(nodes[j]);
				if(p>0 && r>0)
				{
					g(i,j)=g(i,j)+elements[e].Cond()(p,r);
				}else{
					g(i,j)=g(i,j);
				}
			}
		}
	}
	return g;
}

Matrix PMesh::assembleGlobalGradien(const double &kx,const double &ky, const double &w)
{
	Matrix g=Matrix::zeros(numOfN,numOfN);
	for(int i=1;i<=numOfE;i++){
		elements[i].constructGradien(kx,ky,w);
	}
	
	int p=0;
	int r=0;
	
	for(int e=1;e<=numOfE;e++){
		for(int i=1;i<=numOfN;i++)
		{
			p=elements[e].indexOf(nodes[i]);
			for(int j=1;j<=numOfN;j++)
			{
				r=elements[e].indexOf(nodes[j]);
				if(p>0 && r>0)
				{
					g(i,j)=g(i,j)+elements[e].Gradien()(p,r);
				}else{
					g(i,j)=g(i,j);
				}
			}
		}
	}
	return g;	
}

Matrix PMesh::assembleGlobalConv(const double &h, const double &w)
{
	Matrix g=Matrix::zeros(numOfN,numOfN);
	for(int i=1;i<=numOfE;i++){
		elements[i].constructConv(h,w);
	}
	
	int p=0;
	int r=0;
	
	for(int e=1;e<=numOfE;e++){
		for(int i=1;i<=numOfN;i++)
		{
			p=elements[e].indexOf(nodes[i]);
			for(int j=1;j<=numOfN;j++)
			{
				r=elements[e].indexOf(nodes[j]);
				if(p>0 && r>0)
				{
					g(i,j)=g(i,j)+elements[e].Conv()(p,r);
				}else{
					g(i,j)=g(i,j);
				}
			}
		}
	}
	return g;	
}

Matrix PMesh::assembleStiff(const double &rho,const double &cp, const double &kx, const double &ky, const double &h, const double &w){
	
	Matrix global=Matrix(numOfN,numOfN);
	Matrix cond=assembleGlobalCond(rho,cp);
	Matrix conv=assembleGlobalConv(h,w);
	Matrix gradien=assembleGlobalGradien(kx,ky,w);
	
	Matrix k=gradien+conv;	
	global=cond+(k)*(beta*dt);
	return global;
	
}

Vector PMesh::assembleGlobalForce(const double &q,const double &h,const double &w,const double &gen, const double &Tf){
	
	Vector global = Vector::zeros(numOfN);
        for (int i = 1; i <= numOfE; i++) {
		elements[i].constructForce(q,h,w,gen,Tf);
	}

	int p = 0;
	for (int e = 1; e <= numOfE; e++) {
		for (int i = 1; i <= numOfN; i++) {
			p = elements[e].indexOf(nodes[i]);
			if (p > 0) {
				global(i) = global(i) + elements[e].Force()(p);
			} else {
				global(i) = global(i);
			}

		}
	}
	return global;
}
Vector PMesh::assembleForce(const double &q,const double &h,const double &w,const double &gen, const double &Tf,const double &rho, const double &cp, const double &kx, const double &ky){
	
	Vector force=assembleGlobalForce(q,h,w,gen,Tf);
	Vector g=Vector(numOfN);
	Vector T=nodesT();
	Matrix cond=assembleGlobalCond(rho,cp);
	Matrix conv=assembleGlobalConv(h,w);
	Matrix gradien=assembleGlobalGradien(kx,ky,w);
	Matrix k=conv+gradien;
	
	g=(cond-(k*(1-beta)*dt))*T + force*dt;
	return g;
}

Vector & PMesh::nodesT()
{
	Vector *v=new Vector(numOfN);
	for(int i=1;i<=numOfN;i++){
		v->set(i,nodes[i].T());	
	}
	return *v;
}

void PMesh::readNodeData(const char *filename){
	ifstream is(filename,ios::in);
	double x,y;
	int i=1;
	std::cout<<"reading node data...\n";
	while(!is.eof() && i<=numOfN){
		is>>x;
		is.setf(ios::skipws);
		is>>y;
		nodes[i].setLabel(i);
		nodes[i].setX(x);nodes[i].setY(y);
		i++;
	}
	
	for(int i=1;i<=numOfN;i++){
		std::cout<<"Node "<<nodes[i].Label()<<" : "<<nodes[i].X()<<","<<nodes[i].Y()<<"\n";
	}
	
	is.close();
}

void PMesh::readElementData(const char *filename){
	ifstream is(filename, ios::in);
	int i,j,k;
	int e=1;
	std::cout<<"\nreading element data...\n";
	while(!is.eof() && e<=numOfE)
	{
		is>>i;
		is.setf(ios::skipws);
		is>>j;
		is.setf(ios::skipws);
		is>>k;
		is.setf(ios::skipws);
		
		elements[e].setLabel(e);
		elements[e].addNode(1,Node(nodes[i]));
		elements[e].addNode(2,Node(nodes[j]));
		elements[e].addNode(3,Node(nodes[k]));
		
		e++;
	}
	is.close();
	
	for(int i=1;i<=numOfE;i++){
		std::cout << "element " << i << " : " << elements[i][1].Label() << ","
				<< elements[i][2].Label() << ","
				<< elements[i][3].Label() << "\n";
	}
	
}
